Polynomial reconstruction in SharpClaw by sanromd · Pull Request #492 · clawpack/pyclaw

sanromd · 2014-12-28T11:18:44Z

Polynomial reconstructions for SharpClaw (non-dissipative option for EM)

This PR adds linear polynomial reconstructions to SharpClaw. The main motivation is electromagnetic wave problems, where a reconstruction that adds no numerical dissipation is useful. The branch dates from 2014 and has been brought up to date with master (merges cleanly); the reconstruction itself has been rewritten, as described below.

What it adds

solver.lim_type gains three values. All take an even solver.poly_order N in {4, 6, 8, 10}, with a default of 4:

`lim_type`	name	cells per face value	formal order	numerical dissipation
4	`'poly'`	N, centred on the face	N	none
5	`'poly_dissipative'`	N, shifted one cell upwind	N	strongest
6	`'poly_cell_centered'`	N+1, centred on the cell	N+1	moderate

With 'poly', both sides of every face get the same value, so the Riemann solver sees no jump. For linear, constant-coefficient advection, the semi-discrete scheme has Re λ = 0 at every wavenumber, so the only damping comes from the time integrator.
'poly_cell_centered' with N = 4 reproduces WENO5's ideal linear weights. In general it equals 'poly' of order N+2 plus an upwind term proportional to the face jump (an exact identity).
They are Fortran-only, component-wise and unlimited, so they suit smooth solutions and oscillate at discontinuities. They work in 1D, 2D and 3D: flux2/flux3 call flux1 for each slice.
num_ghost is set automatically to N/2 + 1.

API

solver.lim_type accepts the three polynomial names or their integer values: 'poly' (4), 'poly_dissipative' (5) and 'poly_cell_centered' (6). It is still stored as an integer.
solver.lim_type = 'poly' selects fourth-order face-centred reconstruction by default. Set solver.poly_order to choose a different polynomial order.
solver.weno_order and solver.poly_order are independent. Setting either leaves lim_type alone, and switching methods preserves both orders.
An invalid poly_order raises ValueError during solver setup when a polynomial reconstruction is selected. Previously it hit a Fortran stop, which killed the interpreter.

Behaviour changes for existing code

None for existing WENO configurations: fifth-order WENO remains the default, solver.weno_order works without a deprecation warning, and lim_type=3 uses master's ghost-cell count.
The 1D example setup() signatures are unchanged. setup(weno_order=11) still runs WENO11 on lim_type=2 and reproduces the existing regression value.

Documentation and reproducibility

development/poly_reconstruction/ contains:

poly_weights.py: generates src/pyclaw/sharpclaw/poly.f90 from the exact weights (--write, --check, --table).
A sine-function example and toy tests.
A README with figures covering the stencils, the primitive-function construction, convergence, and dissipation/dispersion, including phase-speed and damping tables at 4–16 points per wavelength.

Tests

src/pyclaw/tests/test_sharpclaw_poly.py covers:

reconstruction order for all 12 combinations, and face single-valuedness for 'poly';
solver convergence, energy conservation, and stability at N = 8 and 10;
a 2D run;
the polynomial lim_type names, the default orders, and the independence of weno_order, poly_order and method selection;
a check that poly.f90 matches its generator and that the compiled weights equal the exact fractions to 1e-15.

After the API update, 75 targeted tests passed locally with rebuilt Fortran kernels, including the polynomial tests and the existing 1D acoustics and advection regressions.

🤖 Generated with Claude Code

Historical

The original proposal from 2014 was:

This PR adds a routine for polynomial based reconstruction in SharpClaw. In detail,

4th and 6th order Polynomial reconstrcution

8th and 10th order reconstruction are also included but have not been fully tested yet

renames "weno_order" to "reconstruction_order"

deallocate calls for WENO5 and Poly_comp

sets number of ghost cells for either Weno or Poly in solver

documentation

ketch · 2014-12-30T09:59:10Z

Just a reminder for anybody testing this -- you'll need to rebuild:

cd clawpack
pip install -e .

ketch · 2014-12-30T10:18:31Z

Just a note: renaming weno_order to reconstruction_order is backward-incompatible, but it is a small enough change that I think it's okay for this to go into 5.3.

@ketch, if needed for 5.4 I can add backward-compatibility...

Maybe that's a good idea; we can just emit a warning saying that the old attribute will be deprecated in future versions.

ketch · 2014-12-30T10:24:33Z

I just merged in the fix for the doctests (from master) so this is passing now.

ketch · 2014-12-30T10:27:41Z

This should say

From 5th to 17th order (odd) for WENO reconstruction and from 4th to 10th order (even) for piecewise polynomial reconstruction.

Ah, I just noticed that it's wrong in master...

ketch · 2014-12-30T10:30:15Z

This looks great and will be very useful. It would of course be nice to have:

A test that uses the new reconstruction
At least some comments at the top of poly.f90 documenting the functions there.

coveralls · 2014-12-30T10:46:16Z

Coverage remained the same when pulling 5f6ee25 on sanromd:dev-poly into e20a787 on clawpack:master.

sanromd · 2014-12-30T15:30:26Z

@ketch working on

acoustics_1d test
comments on poly
documentation on solver

mandli · 2015-10-24T14:25:00Z

What's the status of this PR?

sanromd · 2015-10-24T20:32:55Z

I made some modifications to it in my fork and rebased it... I haven't got
time to finalize testing. I will update it soon.

On Saturday, October 24, 2015, Kyle Mandli notifications@github.com wrote:

What's the status of this PR?

—
Reply to this email directly or view it on GitHub
#492 (comment).

Damián

ubi dubium ibi libertas

mandli · 2015-10-24T23:32:42Z

Thanks for the update. Would you like this included in the 5.3.1 release or would this be better for the 5.4.0 release?

ketch · 2015-10-25T14:19:02Z

If @sanromd can add a test or two soon, it would be great to get this into 5.4.0. It probably shouldn't go into 5.3.1 since it's a major feature and 5.4 is coming very soon anyway.

coveralls · 2016-03-18T10:25:02Z

Changes Unknown when pulling 458933d on sanromd:dev-poly into * on clawpack:master*.

coveralls · 2016-10-05T09:56:53Z

coverage: 49.588% (+2.5%) from 47.079% — sanromd:dev-poly into clawpack:master

…order for Poly is >6

…ruction can be selected using lim_type=3

…ruction_orders included'

mandli · 2026-10-13T13:11:25Z

Hi @sanromd — we're working through long-open Clawpack PRs ahead of v5.15. The polynomial reconstruction here is still of interest, but master has moved a long way since 2014, and the reconstruction code in SharpClaw needs some rework of its own (we are planning to look at the WENO/characteristic-decomposition pieces). Are you still interested in this? If so, the path would be to rebase onto current master and add tests for the 4th- and 6th-order reconstructions; the untested 8th/10th orders could be split out. If not, we'll keep the branch and fold the ideas into that rework with credit to you.

sanromd · 2026-10-14T04:58:41Z

Hi @mandli! I would love to support this effort and PyClaw. The polynomial reconstruction is very useful for the time-varying and heterogeneous Maxwell’s equations, which I have also been planning on submitting as a PR. I am at conference this week, I can work on the polynomial reconstruction this weekend following your advise. Sent from Gmail Mobile

…

On Tue, Oct 13, 2026 at 16:11, Kyle Mandli ***@***.***> wrote: *mandli* left a comment (clawpack/pyclaw#492) <#492 (comment)> Hi @sanromd <https://github.com/sanromd> — we're working through long-open Clawpack PRs ahead of v5.15. The polynomial reconstruction here is still of interest, but master has moved a long way since 2014, and the reconstruction code in SharpClaw needs some rework of its own (we are planning to look at the WENO/characteristic-decomposition pieces). Are you still interested in this? If so, the path would be to rebase onto current master and add tests for the 4th- and 6th-order reconstructions; the untested 8th/10th orders could be split out. If not, we'll keep the branch and fold the ideas into that rework with credit to you. — Reply to this email directly, view it on GitHub <#492?email_source=notifications&email_token=AAHAOU2DUTYMQTD5ZX7VCLT5TYSZLA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTMMJUGYZTENJRGE22M4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-6146325115>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAHAOU4YISQ3XNQZBRUNGL35TYSZLAVCNFSNUABBKJSXA33TNF2G64TZHMYTMMRYG4YTCO2JONZXKZJ3GUZDSNZYG4YTBILWAI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

Bring dev-poly up to date with clawpack/pyclaw master (78 commits). Conflict resolution: the five 1D examples/tests take master's version (pytest-style tests, num_cells/nx defaults, RK coefficients) with weno_order renamed to reconstruction_order, which is the only name the dev-poly SharpClaw solver recognises. acoustics_1d keeps dev-poly's lim_type argument. Also add poly.f90 to the SharpClaw Fortran extensions in meson.build, since reconstruct.f90 now uses the poly module and upstream builds with Meson rather than setup.py. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

The previous polynomial weights were point-value Lagrange interpolation weights. PyClaw stores cell averages, so every order converged at 2nd order, and with the default three ghost cells orders 8 and 10 read outside the array and blew up. poly.f90 now builds interface values from cell averages with exact weights, in three variants that all take an even reconstruction_order N in {4, 6, 8, 10}: lim_type 4, 'poly' N cells centred on the face; order N; single-valued at faces, so no numerical dissipation lim_type 5, 'poly_dissipative' same width shifted one cell upwind; order N, dissipative lim_type 6, 'poly_cell_centered' N+1 cells centred on the cell; order N+1, dissipative (WENO with ideal linear weights) solver.py sets num_ghost = N/2 + 1 for these, validates the order and kernel language in Python instead of a Fortran stop, and no longer applies the polynomial ghost-cell formula to lim_type 3, which again behaves as upstream. lim_type also accepts names ('godunov', 'tvd', 'weno', 'weno5', 'poly', 'poly_dissipative', 'poly_cell_centered') and is still stored as an integer. weno_order is kept as an alias of reconstruction_order that issues a FutureWarning; setting it raised TypeError on this branch before. WENO code paths are unchanged. Adds tests for reconstruction order, face single-valuedness, solver convergence, energy conservation, high-order stability, a 2D run, the lim_type names and the weno_order alias. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

solver.weno_order and solver.poly_order are now permanent properties that read back the order of the active reconstruction family (None otherwise) and, when set, select that family: weno_order switches to lim_type 2 unless lim_type is already 2 or 3, and poly_order switches to lim_type 4 unless lim_type is already 4, 5 or 6. reconstruction_order stays the canonical attribute. weno_order no longer warns. SharpClawSolver.set_reconstruction() applies lim_type, reconstruction_order, weno_order and poly_order together and raises ValueError for contradictory combinations. The 1D examples that took an order argument (acoustics_1d, advection_1d, advection_1d_nonunif, variable_coefficient_advection, cubic) accept weno_order= and poly_order= again, so upstream-style calls such as setup(weno_order=11) run WENO11 on lim_type 2 as before. Tests cover the properties, the set_reconstruction rules, the example keywords (including the upstream WENO11 regression value), and check that poly.f90 matches its generator in development/ when present. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

development/poly_reconstruction/ contains: - poly_weights.py: exact (rational) weights for lim_type 4, 5 and 6, their Fourier symbols, and the generator of src/pyclaw/sharpclaw/poly.f90 (--write, --check, --table). - sine_reconstruction.py: reconstructs sin(2 pi x) from N cell averages with every type and order; prints the error and observed order and, with --fortran, compares with the compiled code. - test_poly_reconstruction.py: toy tests (accuracy, polynomial exactness, stability and dissipation, the old point-value weights being second order, generator up to date). - make_figures.py and figures/*.png: figures made with NumPy and matplotlib only, run by hand; output is deterministic. - README.md: what the reconstructions compute, why averages must not be treated as point values, the primitive-function construction, the stencils, weight tables, accuracy, and dissipation/dispersion, including the exact identity that 'poly_cell_centered' of order N is 'poly' of order N+2 plus an upwind term proportional to the jump. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

Setting solver.weno_order or solver.poly_order, or passing weno_order= or poly_order= to an example setup() (through set_reconstruction), issues a FutureWarning saying future releases may deprecate them and pointing to reconstruction_order. The warning is attributed to the caller's line (the assignment or the setup() call). Reading the properties stays silent, and reconstruction_order never warns. Tests check the warnings, their location, that reconstruction_order and lim_type work as setup() keywords in every example without warnings, and move the generator fixture to module scope (pytest deprecates class-scoped fixtures defined as instance methods). The README presents reconstruction_order as the keyword to use. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

Each paragraph and list item is now a single line, so viewers that render single newlines as line breaks no longer break sentences in the middle. The text is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

- weno_order no longer switches a TVD or Godunov solver to WENO; it only moves a polynomial reconstruction back to lim_type 2, as it was inert on non-WENO solvers before. Its getter again reads 5 on those solvers. - setup() validates lim_type, the WENO order (odd, 5-17) and the kernel support for lim_type -1/0 in Python, raising ValueError or NotImplementedError. These cases used to reach a Fortran 'stop', which ends the process with exit status 0 and no exception. - set_reconstruction() takes a stacklevel argument so a direct call attributes its FutureWarning to the caller; the examples pass 4. - poly.f90 (via its generator) uses 'error stop' for its own checks. - examples/euler_1d/Makefile lists poly.f90, which reconstruct.f90 uses. - cubic.py keeps upstream's positional order (order, N); lim_type is appended after N. - Regression tests for each item; 115 tests in test_sharpclaw_poly.py. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

mandli · 2026-10-14T12:34:55Z

Thanks @sanromd, let us know when you have something you want us to take a look at.

Brings in the CI changes on master (clawpack#763, clawpack#768-clawpack#771: macOS moved to a nightly run, PR tests pinned to ubuntu-24.04, Dependabot, action bumps) and the petclaw six removal (clawpack#765). No file overlaps with this branch. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01ThxEJQ7zPNDJ2WipNvNBpT

sanromd · 2026-11-06T10:51:59Z

HI @mandli, greetings.
I teamed up with Claude Code to do the code review and update, given my limited time. We went through the work together and updated the earlier polynomial reconstruction.

The latest update is current with master (2630ce4). CI passes on both ubuntu-24.04 jobs (opt and debug), and GitHub shows no conflicts.
Tests: src/pyclaw/tests/test_sharpclaw_poly.py covers the reconstruction order for every type and order combination. It also checks solver convergence, energy conservation for the non-dissipative option, a 2D run, and the API and warnings. A further test confirms that poly.f90 matches the script that generates it.
weno_order still works and now issues a FutureWarning that points to reconstruction_order, as discussed in 2015. The docstring for the order range has been updated.
Existing behaviour: No WENO routine was changed, and existing WENO regression results are unchanged.

mandli · 2026-11-06T13:30:30Z

Thanks @sanromd! @ketch can you review this?

ketch · 2026-11-08T12:40:10Z

 from clawpack import riemann


 def setup(use_petsc=False, kernel_language='Fortran', solver_type='classic',


I don't like the idea that our example code explicitly names three different arguments (weno_order, poly_order, and reconstruction_order) that are incompatible with each other. Also, it seems to set them all to zero, so this example won't actually run if called with default arguments, right?

I think I would rather just make the examples use the new interface. It's nice to keep backward compatibility for users' code, but the examples should show the right way to do things in the current version of the code.

ketch · 2026-11-08T13:26:37Z

I decided to keep all of the existing interface and simply add the new lim_type options. I've just pushed a patch that does this. @sanromd please take a look and tell me if it seems okay to you.

It might be nice to consolidate and overhaul the reconstruction options API, but that would be part of a major update.

ketch reviewed Dec 30, 2014
View reviewed changes

sanromd force-pushed the dev-poly branch from 5f6ee25 to 48952d2 Compare November 4, 2015 10:35

sanromd force-pushed the dev-poly branch from 48952d2 to 458933d Compare March 18, 2016 10:09

sanromd force-pushed the dev-poly branch from 458933d to c0095d2 Compare October 5, 2016 09:40

sanromd force-pushed the dev-poly branch from c0095d2 to cceb0d7 Compare July 26, 2017 12:32

sanromd added 12 commits October 20, 2020 11:11

polynomial interpolation routines added

1d44160

solved merge conflicts by changing weno_order to interpolation_order

62ecd18

added calls to module poly in resonsctruct and set setup to compile poly

f9edcb1

added WARNING and ERROR messages to reconsctruct when reconstrcution …

ffe513f

…order for Poly is >6

added poly to custom sharpclaw build in euler_1d

9a648ee

solved merge conflict by changin weno_order to interpolation_order

1dbfabe

added allocated, deallocate and call for polyinterpolation and WENO5

5a53304

allocate in reconstruct

943c2f2

added documentation to solver for lim_type=4 (polynomial)

f7f409f

changed weno_order for reconstruction_order in function setup

e51b1d5

added lim_type to function setup

385bc83

declared self.num_ghosts as int

e5ac2fe

sanromd force-pushed the dev-poly branch from cceb0d7 to e5ac2fe Compare October 20, 2020 09:24

sanromd added 4 commits October 20, 2020 14:55

removed case 4 for lim_type select in reconstruct. Polynomial reconst…

4de6e91

…ruction can be selected using lim_type=3

corrected error message in poly routine; it now describes the reconst…

0f939e8

…ruction_orders included'

Merge branch 'clawpack:master' into dev-poly

d572d19

Merge branch 'clawpack:master' into dev-poly

8b7ad7d

claude added 6 commits October 14, 2026 05:17

sanromd force-pushed the dev-poly branch from e011044 to c4d27ec Compare October 14, 2026 08:33

sanromd force-pushed the dev-poly branch from 0728288 to 45d1b51 Compare October 14, 2026 10:56

ketch reviewed Nov 8, 2026

View reviewed changes

Keep WENO and polynomial reconstruction orders independent

73f22d0

		from clawpack import riemann


		def setup(use_petsc=False, kernel_language='Fortran', solver_type='classic',

Conversation

sanromd commented Dec 28, 2014 • edited by ketch Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Polynomial reconstructions for SharpClaw (non-dissipative option for EM)

What it adds

API

Behaviour changes for existing code

Documentation and reproducibility

Tests

Historical

Uh oh!

ketch commented Dec 30, 2014

Uh oh!

ketch Dec 30, 2014

Choose a reason for hiding this comment

Uh oh!

sanromd Nov 4, 2015

Choose a reason for hiding this comment

Uh oh!

ketch Nov 4, 2015

Choose a reason for hiding this comment

Uh oh!

ketch commented Dec 30, 2014

Uh oh!

ketch Dec 30, 2014

Choose a reason for hiding this comment

Uh oh!

ketch Dec 30, 2014

Choose a reason for hiding this comment

Uh oh!

ketch commented Dec 30, 2014

Uh oh!

coveralls commented Dec 30, 2014

Uh oh!

sanromd commented Dec 30, 2014

Uh oh!

mandli commented Oct 24, 2015

Uh oh!

sanromd commented Oct 24, 2015

Uh oh!

mandli commented Oct 24, 2015

Uh oh!

ketch commented Oct 25, 2015

Uh oh!

coveralls commented Mar 18, 2016

Uh oh!

coveralls commented Oct 5, 2016 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mandli commented Oct 13, 2026

Uh oh!

sanromd commented Oct 14, 2026 via email

Uh oh!

mandli commented Oct 14, 2026

Uh oh!

sanromd commented Nov 6, 2026

Uh oh!

mandli commented Nov 6, 2026

Uh oh!

ketch Nov 8, 2026

Choose a reason for hiding this comment

Uh oh!

ketch commented Nov 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

sanromd commented Dec 28, 2014 •

edited by ketch

Loading

coveralls commented Oct 5, 2016 •

edited

Loading

ketch commented Nov 8, 2026 •

edited

Loading